» % * % 


NASA Contractor Report 187441 
ICASE Report 90-78 


ICASE 


THE ROUTE TO CHAOS FOR THE 
KURAMOTO-SIVASHINSKY EQUATION 


Demetrios T. Papageorgiou 
Yiorgos S. Smyrlis 


Contract No. NAS 1-18605 
October 1990 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, Virginia 23665-5225 

Operated by the Universities Space Research Association 


ms/\ 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 5225 

(• „ '-1 . 7 L ) I I. KTjT- Tn CH4U-1 T o- 

K ! j * ft ' i -r-iVA~ntN5AV r ^!ATI r; Hn.ll - Y° r ! 



The route to chaos for the Kuramoto-Sivashinsky equation.' 


Demetrios T. Papageorgiou 
Department of Mathematics and 
Center for Applied Mathematics and Statistics 
New Jersey Institute of Technology 
Newark, New Jersey 07102 

Yiorgos S. Smyrlis 
Department of Mathematics 
University of California at Los Angeles 
Los Angeles, California 90024-1555 


ABSTRACT 

We present the results of extensive numerical experiments of the spatially periodic ini- 
tial value problem for the Kuramoto-Sivashinsky equation. Our concern is with the asymp- 
totic nonlinear dynamics as the dissipation parameter decreases and spatio-temporal chaos sets 
in. To this end the initial condition is taken to be the same for all numerical experiments (a 
single sine wave is used) and the large lime evolution of the system is followed numerically. 

Numerous computations were performed to establish the existence of windows, in parameter 
space, in which the solution has the following characteristics as the viscosity is decreased: A 
steady fully modal attractor to a steady bimodal attractor to another steady fully modal attrac- 
tor to a steady trimodal attractor to a periodic (in time) attractor, to another steady fully 
modal attractor, to another time-periodic attractor, to a steady tetramodal attractor, to another 
time-periodic attractor having a full sequence of period-doublings (in the parameter space) to 
chaos. Numerous solutions are presented which provide conclusive evidence of the period- 
doubling cascades which precede chaos for this infinite-dimensional dynamical system. These 
results permit a computation of the lengths of subwindows which in turn provide an estimate 
for their successive ratios as the cascade develops. A calculation based on the numerical 
results is also presented to show that the period doubling sequences found here for the 
Kuramoto-Sivashinsky equation, are in complete agreement with Feigenbaum’s universal con- 
stant of 4.669201609... . Some preliminary work shows several other windows following the 
first chaotic one including periodic, chaotic and a steady octamodal window; however the 
windows shrink significandy in size to enable concrete quantitative conclusions to be made. 

* Research was supported in part by the National Aeronautics and Space Administration under NASA 
Contract No. NAS 1-18605 while the authors were in residence at the Institute of Computer Applications in Sci- 
ence and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. Additional support for 
the second author was provided by ONR Grant N-00014-86-K-0691 while he was at UCLA. 
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1. Introduction. 

The present study is concerned with the long time behavior of the one-dimensional Kuramoto- 
Sivashinsky (KS) equation, conveniently normalized as 

u, + uu x + u„ + = 0 , 

(x,t) e R l xR + , (1.1) 

w(x.O) = u 0 (x) , u(x+2n,t) - u(x,t) . 

The parameter v is positive and takes the role of the viscosity of the system. 

The long time behavior of the KS equation is characterized by the negative - therefore destabiliz- 
ing - second order diffusion, the positive - therefore stabilizing - fourth order dissipation and the non- 
linear coupling term. As soon as the value of v decreases below 1 at least one Fourier frequency 
becomes unstable in the linearized model. Nevertheless the solution does not grow exponentially fast 
as the linear analysis around the null state would suggest but instead stays bounded due to the transfer 
of energy from the linearly unstable low modes, the number of which is mod(v _l/2 ), to the linearly 
stable high modes. The previous statement has been observed numerically and also established analyti- 
cally (see Constantin et al. [8].) 

A very attractive feature of (1.1), is that even though it appears to be a very simple partial 
differential equation it produces complicated dynamics in both space and time (chaotic behavior). The 
equation arises in a variety of physical problems, such as in chemical physics for propagation of con- 
centration waves (Kuramoto [24], Kuramoto and Tsuzuki [25], [26]), plasma physics (Cohen et al. [6]), 
flame propagation (Slvashinsky [35]), reaction diffusion combustion dynamics (Sivashlnsky [36]), free 
surface film-nows (Benney [4], Sivashinsky and Michelson [37], Shlang and Sivashinsky [33], Hooper 
and Grimshaw [17]), and two-phase flows in cylindrical geometries (Frenkel et al. [13], Papageorgiou 
et al. [31]). 

There have been numerous computational studies of (1.1) (Cohen et al. [6], Sivashinsky and 
Michelson [37], Aimar [1], Manneville [29], Frisch, She and Thual [14], Hyman and Nicolaenko [18], 
Hyman, Nicolaenko and Zaleski [19], Kevrekldis, Nicolaenko and Scovel [22], Greene and Kim [15] 
to mention a few), as well as analytical ones ( Constantin et al. [8], Foias et al. [12], Babin and Vishik 
[3], etc.). The overall picture that emerges is that as the viscosity parameter v is decreased, the large 
time behavior of the system varies from steady-stale (for example cellular) solutions to complicated 
time-oscillatory (chaotic) ones. Such behavior has led to the conjecture that the solution behaves 
according to low modal dynamics with the determining number of Fourier modes being proportional to 
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v- ,/2 (i.e. Pomeau and Manneville [32]). Significant analytical results have provided estimates for the 
dimension and the nature of the global attractor (see [8],[12]), but the number of determining modes 
has not yet been estimated mathematically with satisfactory accuracy and thus we chose initially to 
integrate the partial differential equation directly rather than work with a truncated system. A truncated 
system of coupled nonlinear ODE's for the Fourier coefficients was also integrated and the number of 
determining coefficients was established experimentally. Since the dynamic behavior of the system is 
of a low modal structure (see earlier references), we found that all the Important phenomena, for exam- 
ple period doubling cascades to chaos, are completely captured by a truncated system - this was 
verified by comparison of numerical solutions of the infinite dimensional dynamical system (PDE) and 
the truncated finite dimensional one, with emphasis on a window that yields a sequence of eight period 
doublings before chaos sets in. From hindsight we discovered that in the cases of the smaller values 
of v studied here, as many as approximately 40 coupled nonlinear equations must be retained if the 
time-periodic solutions are to be described with satisfactory accuracy. 

Of particular relevance to the present work are the studies of Hyman and Nicolaenko [18] (also 
Hyman et al. [19], Kevrekides et al. [22], and Jolly et al. [20]). These studies are mostly concerned 
with the extensive numerical study of the related KS equation for the integral in space of our function 
u ( x ,t ). The equation used is 

«/ + 4« ira + a(«„ + ju/) = 0 , 

(x,t) e R*xR + , (1.2) 

u(x,0) = u 0 (x) , u(x+2n.t) = u(x,t) . 

The parameter a above is related to v by a = 4/v. The range of a used was 4 £ a £ 120, which 
corresponds to 1/30 £ v £ 1. Various attractors were constructed by changing the initial condition as 
well as a, by continuation methods where the initial condition is renormalized to its current value 
before a is changed. The bifurcation diagram reported contains windows of various attractors (fully 
modal, bimodal, trimodal and tetramodal fixed points) separated by oscillatory or chaotic orbits. A 
time-periodic orbit is reported in a small window 17.30 £ a £ 22.50. 

The current study obtains a bifurcation diagram which is analogous to the above, but some essen- 
tial differences must be pointed out. First of all, our concern is not with the full description of the 
inertial manifolds which are by no means easy to obtain since the Initial conditions belong to an 
infinite-dimensional vector-space. Instead we want to investigate transition to chaos, which, as we 
have observed, occurs according to the classical routes of period-doubling bifurcations. Such 
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transi lions to chaos have been described for iteration maps (Feigenbaum [10], Collet and Eckman [7], 
see Devaney [9] for a full description) as well as for low dimensional systems of ordinary differential 
equations such as the forced Duffing equation which is a second order system, and for the Lorenz 
attractor (see Guckenheimer and Holmes [16] and references therein). Numerical evidence of a period- 
doubling sequence has been presented by Knoblock et al. [23] and Lennie et al. [27] for the two- 
dimensional convection in a Boussinesq fluid. Many interesting results have also been found for the 
Ginzburg-Landau equation by Keefe [21] (see also Moon, Huerre and Redekopp [30], Sirovich and 
Rodriguez [35]) including a period-doubling sequence that supports an accurate estimate for the 
Feigenbaum number. Finally evidence of a period-doubling sequence has been discovered for a trun- 
cated model of the Kuramoto-Sivashinsky equation by Brown et al. [5] (see also Jolly et. al [20]); two 
period-doublings are reported in a window which corresponds to the tiny range 
0.121223 £ v £ 0.121304 for our problem. In the results we present here the initial condition is taken 
to be the same for all v. Global attractors are insensitive to such a choice, but in regimes where the 
solution becomes time-periodic and period-doubling bifurcations take place the present choice of initial 
condition enables both an accurate large-time integration as well as a quantification of the cascades 
(see later also). We emphasize here that in this article we exclude a complete study of the phase-space 
of initial conditions. 

The structure of the paper is as follows. Section 2 describes the numerical methods used, while 
in Section 3 we begin to describe each window separately by showing representative numerical experi- 
ments in each case. Here we describe the first five windows which support solutions that are steady (i) 
(fuUy modal, bimodal, fully modat again, trimodal), time-periodic (ii), fully modal steady (iii), time- 
periodic again (iv) and steady tetramodal (v). In Section 4 results for the sixth window and beyond 
are given; the sixth window is time-periodic and contains a complete infinite sequence of period- 
doubling bifurcations leading to chaos. Results are presented for both a direct simulation (classical 
many mode discretization) and a truncated model (traditional Galerkin) with the low modal behavior of 
the dynamical system clearly brought out. The results are used to compute the Feigenbaum number for 
the cascade which is found to be in excellent agreement with Feigenbaum ’s universal theory. This 
Section also includes some preliminary results of solutions at smaller v showing the appearance of 
steady attractors of decreasing spatial period; time-periodicity and chaotic behavior persist as expected. 
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2. Numerical methods. 

As mentioned earlier, our concern here is primarily with the numerical simulation of (1.1). Two 
numerical methodologies are used: A classical numerical discretization with many degrees of freedom 
(which we term the "direct" method) and a traditional Galerfcin projection onto a low-dimensional iner- 
tial manifold that yields a coupled nonlinear system of ODE’s to solve. In the former case we choose 
to perform the time integration by treating the linear part of the operator in Fourier space, and the non- 
linear term in real space by a scheme that does not add numerical viscosity. This split step method 
maintains second order accuracy by the repeated application of the semi-groups, as suggested by the 
Strang split method. The time-step is chosen adaptively subject to several controls that ensure both 
numerical stability and high accuracy; the error in the experiments reported here is kept well below the 
10 -6 level. Such requirement is necessary to capture the nonlinear coupling of Fourier modes that are 
undergoing period-doubling cascades. For the Galerkin projection the Fourier coefficients were com- 
puted for a pre determined dimension; a fourth order Runge-Kutta method was used to advance the 
nonlinear part in time, while the linear part was treated exactly by use of the split method described 
above together with repeated semi -group application that maintains overall accuracy. A natural accu- 
racy test, throughout the computation, is the conservation of the spatial integral of u . The initial condi- 
tion used is 

u 0 (x) = -sin(x) , 

and a consequence of this is that the integral over a spatial period vanishes for all time, and further 
more the solution maintains its odd parity. More generally, if the spatial profile of the solution is ini- 
tially of the form 

u (x ,0) = £ a k sin kx , 

*21 

then the mechanisms of the evolution of the partial differential equation respect this structure since, 

«« + v iw = - £ (* 2 - v* 4 ) a k sin kx 
*21 

and the nonlinear part becomes, 

uu x = a k sin kxf = 

M *21 

= 1 • • = - |(5>*+xax “ 1 £ fl X a n)) sin **. 

*21 X2l X+Jb=* 

Thus the spatial profile remains a linear combination of sines for all times, l.e. 
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u(x,t) = a k (t) sin kx . 

*2 1 

Due to the above observation the PDE (1.1) becomes equivalent to an infinite dimensional system of 
ordinary differential equations for a k 's , i.e. 

a k = ( k (i) 2 * - vk 4 ) a k + “ j Z a \ a \i) forkZl. (2.1) 

The odd-parity property was utilized by a suitable filtering of Fourier coefficients; the filtering toler- 
ance was set to a value between 10 -10 and 1(T 12 . In the case of the direct computation the spatial 
discretizations were done on grids ranging from 128 to 1024 points. For windows other than the first 
steady attractor, at least 512 points were used. Results that are obtained by the projection method will 
be pointed out together with the number of modes retained for the dynamical system. 

The results presented here were produced after several hundreds of numerical experiments which 
tracked the time evolution of the profile. Us Fourier modes, its energy (L 2 -norm) and the phase-plane 
of the energy and in some marginal cases of certain Fourier coefficients. The time evolution of the 
energy is very representative of the structure of the various attractors. Steady-state attractors evolve to 
a constant large-time state, while periodic or chaotic solutions have a time-oscillatory energy profile. 
Near window and subwindow (period-doubling) boundaries for the smaller values of v at least, the 
time integration required as many as 200 time units (8xl0 5 time steps) in order to sharply classify the 
v-windows and subwindows. AU periodic solutions reported here were derived by a phase-plane 
analysis of the data. As is described in the next Section many period-doubling bifurcations are evi- 
denced from the phase-plane, which Is used in conjunction with energy plots to build the geometrical 
structure of the solution in one time-period and so accurately estimate the the v-position of the hyper- 
bolic point corresponding to a period-doubling bifurcation, for example. Table 1 presents a summary 
of all windows discovered so far, together with sharp computed estimates of window boundaries. 

3. Detailed numerical experiments. Pre-chaotlc windows. 

When the parameter v decreases below unity, the solution undergoes its first Hopf bifurcation. In 
this Section we give a description of the dynamics in each of the first five windows. 

(i) Steady attractors , .05985 £ v < 1 . 

The main characteristic of this window, is the attraction of the initial profile to a steady-state. The 

steady-state is approached exponentially fast and is numerically achieved after a relatively short time 
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(5-10 time units in most cases). However the numerical convergence to the steady steady state is 
significantly delayed towards the lower end of the window where the first periodic attractor sets in. We 
subdivide the window into the following four subwindows : 

(a) A fully modal steady attractor when .248 <> v < 1 . 

(b) A bimodal steady attractor (i.e. where all the odd frequencies vanish) when .0758 <, v £ .247 . 

(c) Another fully modal steady attractor when .06697 £ v £ .0758 . 

(d) A trimodal steady attractor (where only the frequencies which are multiples of three survive) when 
.05985 <.v<, .06695 . 

For values of v near the lower boundary in fact, we find that there are about 30 modes which 
have an amplitude larger than KT 6 . In Figures (1.1), (1.2) and (1.3) we give sets of these attracting 
profiles for different decreasing values of v; Figures (1.1a) and (1.1b) cover ranges of v in the first 
subwindow. Figure (1.2a) covers a range of v at the lower end of the first subwindow as well as at 
the upper end of the second subwindow indicating a smooth transition from the fully modal first 
subwindow to the bimodal second subwindow. Figure (1.2b) also covers a range of v in the second 
subwindow, whereas figure (1.3a) covers a range in the third subwindow where the attracting profiles 
are trimodal. Comparison of figures (1.2b) and (1.3a) provides the indication that the transition of the 
attracting profiles between the two subwindows is smooth. Finally figure (1.3b) covers a range of 
values of v in the fourth and last subwindow where the attracting steady profiles are trimodal. The 
change in the shape of the limiting profiles in this subwindow is hardly visible, even though v changes 
"significantly" from .06695 to .06 . This is not surprising since under a similarity rescaling (see below) 
they correspond to fully modal steady profiles for v between .54 and .603 where not much change is 
observed either. 

For values of v which are near 1, the profile is qualitatively similar to a sine wave. As v 
decreases, however, one more maximum and one more minimum are obtained after the profile attains 
an inflection at the anti symmetry line x = n. Such solutions were also obtained by Frisch et al. [14] 
for values of v larger than or equal to 0.35; as a result of this they did not see profiles with four turn- 
ing points. It is interesting at this point to describe an important property of the steady KS equation; if 
u(x) is a solution corresponding to a viscosity v of the ODE obtained by dropping the time depen- 
dence in (1.1), then for any integer N there exists another solution Nu(Nx) with viscosity v/N 2 . Due to 
the periodic boundary conditions used here, N must be an integer. This self-similar property of the 
ordinary differential equation, satisfied by the steady states of the KS equation, permits the generation 
of a one parameter family of steady-states which grow in amplitude and decrease in period. In sub- 
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section (v), the tetramodal attractor, we give a numerical example of this self-similar property connect- 
ing steady-states for N = 4. 

(ii) Time periodic attractors. .055238 £ v £ .0599 . 

The large time evolution of our Initial value problem for v lying In the range covered by this 
window, is attraction to a time-periodic orbit. In all the results presented here, the period is defined for 
the energy ( L 2 -noim ) of the solution. It is clear that the solution itself can have half the period of 
the energy but such a case was not observed when the Fourier coefficients were considered individu- 
ally. Periodicity is determined by examination of the energy phase -plane ( E versus dEldt ), and the 
period is evaluated from the evolution of the energy with time. All the solutions in this window are 
fully modal, and are characterized by a time-periodicity (the time period will be denoted by x ) which 
begins at a value of approximately x = 1.0183, increases slowly to x = 1.18 at a value of v = .059 and 
then doubles to x = 2.374 at a value of v = .0559. The doubling is Instantaneous, in parameter space, 
and comes about by a slight disalignment of the maxima and minima of the energy-time graph. Some 
representative results from this window are given in Figures (2.1a,b)-(2.3a,b). The solution undergoes a 
period doubling as is evidenced by comparison of the phase-planes at v = .057 and v = .0559 - see 
Figures (2.1a), (2.2a). Figures (2.3a, b) correspond to v = .0555 and it can be seen that the 4-cycle pro- 
perty of the phase diagram is much more pronounced than at the slightly higher value of v = .0559; 
such drastic changes in the dynamics over small domains of the parameter space are typical and are 
encountered later on also. A summary of the results from this window is contained in Table 2. 

(iii) Steady fully modal attractors. .03965 ^ v £ .055235 . 

Here a new steady window is found with fully modal solutions which do not correspond to a res- 
caled value of the parameter v under the self-similar scaling described earlier. The new element here 
is the extremely slow convergence, still exponential nevertheless, towards the fixed point attractor. At 
such values of v the equations must be integrated to long times; this is to be expected since a straight- 
forward perturbation analysis above the window-boundary yields a rate of decay of the distance (for 
example the L 2 -norm) from the attractor, which is inversely proportional to the perturbation of v from 
the boundary. This is illustrated in Figures (3.1a)-(3.1c) for v = .0397. Figure (3.1a) shows the varia- 
tion of tlie energy with time, while Figures (3.1b) and (3.1c) show the variation with time of the L°° 
norm and the L°° norm of the gradient of the solution respectively. It can be seen, that initially the 
solution passes through an oscillatory and almost periodic stage, and at a time of approximately 20 
units Its amplitude follows an envelope that takes the profile into the stationary attractor. In fact 
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during the oscillatory stage, the phase-plane of the energy is similar to that of the periodic window 
described next; once the steady-state begins its attraction the phase curve converges to a sink of course. 

(iv) Time periodic attractors. .03735 £ v £ .0396 . 

Periodicity is again determined by examination of the energy phase-plane, and likewise the period 
Is computed from the energy versus time plots. In this window the period begins at a value of approx- 
imately x=.515 when v = .0396, with a very slow convergence to the periodic attractor, and it 
increases smoothly, in parameter space, until it reaches the value x = .684 when v = .03798 with a 
much faster convergence to the the limit circle; the period then suddenly doubles to x = 1.37 at a value 
of v = .03795. As v is decreased further, the period undergoes a new slow increase to reach a value 
x = 1.43 at v = .0377; a small decrease of v to a value .0376 gives a new period doubling with 
x = 2.87. The period again increases slowly to a value x = 3.0 at v = .03736. At this point a slight 
decrease of v to the value v = .037355 gives a period-halving with the energy having a period x = 1.5. 
The last value of v computed in this window is v = .03735, for which the period is x = 1.5. The next 
numerical experiment with v = .037345 is outside this time-periodic window and yields a tetramodal 
steady-state (see later). These results are summarized in Table 3. 

The results of Table 3 were obtained after numerous computations together with a careful exami- 
nation of the oscillatory dynamics by means of phase-plane analysis, and representative solutions are 
given here. Some cases are described in full including the variation and phase-planes of certain Fourier 
modes, while other cases are described by the energy time series amid the energy phase-plane. These 
two plots are enough to conclude periodicity of solutions and their periods. In Figures (4.1a-d) we 
present results for v = .038 which is a representative parameter at the beginning of the window. These 
Figures show the variation of the energy with time, the energy phase-plane, and the phase planes of the 
first and tenth Fourier modes respectively; time periodicity of solutions is clear and we note that the 
Energy phase plane (Figure (4.1b)) is traced out 21 times since the time interval it depicts contains 21 
periods. The time period of Fourier coefficients is the same as that for the energy, a fact that need not 
hold in general (see earlier comments); the reason it holds in the cases presented by this article, is that 
the actual solution u(x,t) has a non-zero mean part in general. 

As has been mentioned already (Table 3), a period-doubling takes place at v = .03795. This is 
evidenced very clearly in the phase-plane of the energy for this case which is given by Figure (4.2a). 
The phase-plane consists of two loops which are initially almost coincident (in fact they were coin- 
cident at v = .038) and this change produces the period-doubling. The corresponding energy variation 
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wlth time can be found in Figure (4.2b), from which the period is calculated to have a value of 1.37. 
The mechanism of the period doubling can also be seen from (4.2b), where every other wave crest and 
trough are reduced/increased slightly thus producing a phase-plane of index two and a period-doubling. 
We note also that there is no abrupt change in the wave-form that undergoes a period-doubling; only 
the value of the period doubles abruptly. Next we reduce the value of v by 5x1 CT 5 to v = .0379. We 
are still in the same subwindow, but the phase plane diagram is now much more pronounced as Figure 
(4.3a) shows. The results for the last v from this subwindow is presented in Figures (4.3b), (4.3c) and 
(4.3d) which correspond to the phase-planes of the energy, the first and sixth Fourier coefficients 
respectively, at a value of v = .0377. 

A reduction of v to the value .0376 brings us into a new subwindow with the solution undergo- 
ing a period doubling. Evidence for this can be found in Figure (4.4a) which depicts the phase-plane 
of the Energy at v = .0376; the time period is x = 2.87; Figures (4.4b) and (4.4c) are also given to 
show the evolution of the sixth Fourier coefficient with time and its corresponding phase-plane. The 
time-periodicity is again clear and similar results were found for all leading order (i.e. large enough for 
graphical purposes) Fourier coefficients. Such results are of course expected but they serve as a strict 
check on the numerical accuracy of the solutions so that we can confidently state that the time oscilla- 
tions are by no means spurious, but in fact persist in a regular and accurate manner. The last case 
computed in this subwindow has v = .03736. The energy evolution and phase plane curves are given 
in Figures (4.5a) and (4.5b). The solution is just about to undergo a period halving and the index two 
property of the phase plane is only just discernible. A reduction of v to the value .037355 gave a 
period halving, and Figure (4.5c) shows the energy phase-plane for the last member of this window 
that was computed, v = .03735. . As expected, the mechanism that produced period-doublings is now 
working the opposite way to produce period-halvings (similar windows were also found by Knoblock 
et al. [23] for a convection problem). As mentioned earlier, a numerical experiment with v = .037345 
produced a tetramodal steady-stale which is representative of the window described next. 

It is also interesting to observe the spatio-temporal evolution of the profiles, and we present this 
for two representative values of v from this window, v = .039 and .0375 respectively. The way this is 
done is by plotting on the same diagram the solution at successive increasing times. To make the evo- 
lution clearer the profiles are shifted vertically at each time interval producing pictures that capture the 
spatio-temporal evolution and in particular the time-periodicity of the solutions. Figure (4.6) 
corresponds to v = .039 while Figure (4.7) shows the results for v = .0375. 
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(v) Steady tetramodal attractors. .037348 £ v £ .0344 . 

We present only one case (v = .0373) from this window since the results are qualitatively similar, 
the difference being in the final amplitude of the solution. For v = .0373, therefore, a steady-state Is 
reached after approximately seven time units. The evolution of the energy is shown in Figure (5.1a) 
from which it is seen that after a few oscillations in time, the solution gets attracted to a tetramodal 
steady-state; this was determined from inspection of the Fourier coefficients. The corresponding profile 
at t — 10 is depicted in Figure (5.1b) which clearly shows the tetramodal character of the solution. 
This means that the solution is spatially periodic and of period jc/2. Further more, this solution is 
obtainable from the self-similar re-scaling (see earlier) of the bimodal steady state at v = .1492 with 
the value of N = 2. This was checked and agreement is excellent. 

4. Complete period-doubling cascade to chaos and computation of the Feigenbaum number. 

(vi) Time periodic attractors. Complete period-doubling sequence. .0299756 g v £ .0343 . 

This window is to our view, the most important found in this work, and for this reason we con- 
centrate on the results it produces rather than try to sharpen and classify the bifurcation diagram at 
higher values of the viscosity (the classification of the bifurcation diagram at high values of v has been 
the subject of extensive numerical and analytical studies by various authors - [18], [19], (20), [21]). The 
importance of this window is two-fold. First, it supports solutions that start as laminar steady-states 
(tetramodal waves) and then become time-periodic (this is also the case for window (iv) above). 
Second, the time-periods of the solution undergo a complete sequence of period-doubling bifurcations 
as the viscosity v is decreased, resulting in an accumulation point beyond which the solution is chaotic. 
Remarkably enough as many as 5 period-doubling bifurcations were accurately computed by the direct 
method and 8 period-doublings for the Galerkin procedure. This enables, for the first time, the accu- 
rate computation of the Feigenbaum number for the Kuramoto-Sivashinsky equation. Other investiga- 
tions of PDE systems found results that were consistent with the universal Feigenbaum constant but no 
conclusive evidence was shown (see for example Knoblock et al. [22] for a convection problem, and 
Jolly et al. [20] for the KS equation). Significant results were obtained by Keefe [21], however, for 
the Ginzburg-Landau equation. After extensive numerical computations of even periodic solutions 
Keefe found a window containing a complete sequence of period-doubling bifurcations; the first five 
bifurcations were accurately obtained and based on the results the Feigenbaum number for the cascade 
was estimated to be approximately 4.4 to 4.7. The numerical method used in [21] was a pseudo- 
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speclral one with a 32 mode discretization. We mention here that these bifurcations take place on very 
small scales, in parameter space, and require a lot of care numerically, as well as an extensive compu- 
tational effort towards a fairly sharp estimation of subwindow lengths (a subwindow should be under- 
stood to be a region in parameter space that separates period-doubling parameter values). We begin this 
section by presenting the results of the direct method by use of at least 512 Fourier modes. Most of the 
evidence from this approach is graphical and serves to show numerous period-doublings through 
phase-plane analysis. We also re-computed this same window by means of a truncated model that 
retains the leading 36 Fourier modes. The results of the two methods agree to at least 4 decimal 
places (by results here we mean the positions in v-space of the various subwindows); a quantitative 
analysis of the results obtained from the 36-mode Galerkin expansion enables the accurate computation 
of the universal Feigenbaum number of the cascade (see later for details). 

The direct computation of the PDE provides results for the first six subwindows over which the 
solution attains five period-doublings. Improvement to a sixth period-doubling would require a substan- 
tial amount of computational resources and is wise at this point to turn to finite-dimensional versions 
of the equations. For example computations that retain 36 Fourier modes are in general at least 8 times 
faster than direct computations; in the absence of sharp. analytical results direct computations are essen- 
tial in gauging the accuracy or sufficiency of truncated models. As mentioned already, 36 modes are 
capable of capturing a complete set of period-doubling bifurcations that precede chaos, and more 
importantly are in complete agreement with the Feigenbaum universal constant. To give an example of 
the numerical restrictions that are present in either direct or finite-dimensional approaches, our direct 
numerical solutions for v = .0299756 gave a time periodic solution with x = 7.04 while at v = .02997 
the solutions appear to be chaotic. In a distance of less than 5.6xl0~ 5 , therefore, numerous period- 
doubling bifurcations take place (to conjecture an infinite number would be too hasty because the 
chaotic solution at v - .02997 can be thought of as an oscillatory solution with a very large period; the 
period Is so large that discourages any attempt at its numerical computation). It can be seen, therefore, 
that we would very soon run out of computational accuracy, since subwindow lengths would be of the 
same order as computer round-off error. We decided, therefore, to concentrate on the first five subwin- 
dows (eight for the 36 mode model) and compute the Feigenbaum number from these. 

For brevity of the description we will restrict our presentation to phase-plane curves of the 
energy. Each case presented was computed with 512 Fourier modes and approximately 4xl0 4 time 
steps (a run takes about 20 minutes CPU time on a CRAY-11 and about 4 hours CPU time on a 
SPARC work-station; all runs use double precision). Accuracy tests were also made for many runs by 
use of 1024 Fourier modes; all results reported here are insensitive to such refinements as the truncated 
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system also shows (see later). We give representative results from each subwindow and show the 
behavior near subwindow boundaries. The steady tetramodal solution gets attracted to a time-periodic 
orbit when v is reduced below a value of .0345. The phase-plane of the L 2 -norm of the solution is a 
1-cycle initially. In Figure (6.1a) the phase plane for v = .0342 is shown over a time interval from 1.0 
to 34.5. It is clear that the solution gets attracted to a time-periodic orbit of index 1; the period is 0.37. 
As v is decreased the period increases smoothly and the 1-cycle topology of the phase-plane persists 
until the first period doubling suddenly appears at a value of v just above .0304. The phase-plane 
corresponding to v = .0304 is not shown because the period-doubling bifurcation has just taken place 
and within graphical accuracy the phase-plane appears to have the same topology as in the previous 
subwindow. Careful inspection of the left side of the curve, however, shows that there are two distinct 
branches; the period of the solution is 0.89. The two branches begin to move apart as v decreases 
further. An example of this change in the phase-plane, which in turn corresponds to a change in the 
oscillatory structure of the solution, is seen in Figures (6.1b-c); the value of v is .0303 and the period 
is the same as before, x = 0.89. The last picture we present from this subwindow is Figure (6. Id) for a 
value v = .03007. The characteristics of the time oscillations are qualitatively the same as before with 
the two loops well separated and with time-period x = 0.88. 

We next decrease the value of v slightly. A new subwindow is found which develops after a 
period-doubling bifurcation at a value of v between .03007 and .030065. The phase-plane curve now 
becomes a 4-cycle and in one complete time oscillation the energy has eight turning points, four of 
which are maxima and four minima. A characteristic phase-plane curve for bifurcation parameters near 
the upper edge of the second subwindow, is Figure (6.2a) which has v = .03006; the time-period is 
t = 1.76. The index 4 property Is evident and appears by a slight disalignment of the phase plane 
branches. We also note that at this value of v the maxima of the energy are noticeably different 
whereas the minima are almost indistinguishable. Other cases from this subwindow are also given in 
Figures (6.2b-d). The value of the viscosity is v = .030015 , .03 , .029998 respectively. The period is 
1.76 and the 4-cycle property, characteristic of this subwindow, is veiy clear. Figure (6.2d), which has 
v = .029998 is at the lower edge of the second subwindow and is representative of the dynamics just 
before a second period-doubling bifurcation. 

Next the value of v was decreased to .0299975; at first sight this v seems to belong to the second 
subwindow but a careftil inspection of one of the loops shows that a period-doubling takes place some- 
where between v = .029998 and .0299975 leading to an estimated time-period of x = 3.52 and phase- 
plane characteristics for this subwindow of index 8. For smaller values of v the phase-plane loops 
begin to drift apart. Figures (6.3a-c) show a sequence of energy phase-plane diagrams for three 
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decreaslng values of v that belong to this third subwindow; these values are 
v = .029995 , .02999 , .0299815. The last computation in this subwindow is at v = .0299815 and the 
period is x = 3.52. 

The next numerical experiment has v = .029981. Once more, carehil inspection of the phase- 
plane shows that a 16-cycle is obtained and hence another period-doubling bifurcation takes place 
somewhere between v = .0299815 and .029981. To set iliis one must concentrate, for example, on that 
region of the phase plane described by 17.5 < E < 18 and dE/dt near zero. Before the period-doubling 
bifurcation takes place there is a single branch traced out in the phase-plane while now two branches 
are distinguishable within the limits of graphical accuracy. At smaller values of v the branches move 
apart and the 16-cycle property Is very clearly seen. Figures (6.4a-c) give representative phase-planes 
corresponding to a sequence of viscosities which span this fifth subwindow. The bifurcation parameters 
are v = .02998 , .029978 , .0299756 respectively. The value v = .0299756 is also the last computation 
from this subwindow; the period is z = 7.04. 

In Figures (6.5a-d) we present results after another period-doubling occurs. The first case, Figure 
(6.5a), has v = .029975 and the phase-plane appears to be a 32-cycle with a period of at least 14.0. 
The complexity of the dynamics becomes more enhanced as v is decreased even further, and it 
becomes very difficult and computationally expensive to try and follow more period-doubling bifurca- 
tions than the ones described here. Figures (6.5b-d) give a flavor of the complicated dynamics that one 
has to resolve, for the values v = .02997 , .02995 , .0299 respectively. 

Computation of the Feigenbaum number 

As mentioned earlier, the results we have from the direct numerical simulation are probably not 
extensive enough for an accurate determination of the Feigenbaum number. The lacking ingredient is 
the sharp estimation of subwindow boundaries which take on a central role in such a computation. As 
a result of this, a straightforward calculation of the ratio of successive subwindow lengths suffers from 
certain inaccuracies. For completeness, however, v c give the results and stress that they can be 
improved by resort to a much more extensive numerical search. If we label the lengths of successive 
subwindows by /, , /= 1 ,<» ihen the following ratios are obtained from the results of Table 4 : 
li/li = 5.42, / 3 // 2 = 11.35 and ( 4/(3 = 2.77. The cascade is still in its early stages but the distance, in 
v-space, of the fourth subwindow lower boundary from the point where the strange attractor seems to 
commence, is of the order of 1CT 6 . Due to the non-sharp classification of the subwindows, the ratios 
are only of the same order of magnitude as the universal Feigenbaum constant. Our data is much more 
extensive for the 36-mode Galerkin approach and they yield excellent results as we explain next. 
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The period-doubling route to chaos is well documented for certain low dimensional models (see 
earlier references), but it is believed that the Feigenbaum constant, originally found for iteration maps, 
is universal and appears in oilier high-dimensional systems also. In what follows we present numerical 
evidence of tills universality for the Kuramoto-Slvashinsky equation. Suppose that period-doublings 
take place at some successive values, v„ say, of the viscosity parameter. The theory predicts an accu- 
mulation of such bifurcation points as n -> and in fact the universal constant is given by : 


lim 

n — > « 



4.6692016... . 


The obvious way to lest our results against the theory is to form the successive ratios indicated above, 
for our computed subwindow boundaries. The results of this approach are given in Table 5. The first 
column of the Table gives our estimate for subwindow boundaries (for example the first two numbers 
in the first column show the beginning and end respectively of the first periodic subwindow; the 
second and third numbers indicate the beginning and end of the second periodic subwindow etc ), the 
second column shows the subwindow length and the third column contains the ratio of successive 
subwindow lengths. Note that the cascade develops as we move down the Table and by the time the 
third period-doubling bifurcation is attained the geometric ratio of successive bifurcation points is 4.27, 
a number already close to Feigenbaum’s universal constant. It is worth noting that subwindow lengths 
are shrinking significantly and the last subwindow that was computed accurately had a length of 
1.3xlO -7 . It is clear, therefore, that the difficulty of the numerical experiments increases as the cascade 
develops and some special continuation algorithms were used towards the expediency of our estimates. 
Such accuracy measures are essential In the production of four successive bifurcation ratios from the 
fourth to the eighth subwindow, which are at an error of at most 4% from the universal constant of 
4.669201609... . Various other tests can also be performed using the data. For example we can pick a 
subwindow which begins at v„ and which has length Then by use of the computed accumulation 
point v c = .029969275 and by assumption of an infinite cascade, it is easy to show that there exists a 
geometric ratio 



which takes the point v„ to the accumulation point v c . The number 8 above, provides then the 
"Feigenbaum number" for the chosen subwindow. The theory tells us that as the value of v„ is chosen 
to be progressively closer to the accumulation point, 8 then converges to the universal constant. This 
version of the cascade theory was also applied to our numerical data. The results are similar to those in 
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Table 5 but they are slightly more accurate; the gain in accuracy is due to elimination of combination 
of errors that a successive ratio inevitably contains. Finally, we also computed the Feigenbaum number, 
based on subwindows which are not successive but have a factor of 4 difference in their time-periods. 
A ratio is again obtained whose square root should converge to the universal constant after an infinite 
number of bifurcations. After eight bifurcation (even less in fact) the numbers obtained are again very 
similar to those of Table 5 and we omit further details here. 

Finally in this Section we present a brief summary of some of the results obtained for values of v 
which take the dynamics beyond the complete sequence of period-doubling bifurcations just described. 
The richness of the attractors is highly enhanced and the regions in parameter-space over which 
changes take place become smaller. In this article we choose to describe three additional windows; 
these are the chaotic window entered immediately after the period-doublings described above, the 
second is a time-periodic attractor, and the third is a region of parameter space that supports chaotic, 
steady and time-periodic solutions. 

(vii) Window of Chaotic Trajectories. .0252 ^ v £ .029755 . 

No recognizable patterns are observed either in the energy time series or in the phase planes for 
values of the parameter in this window. It should be noted, however, that all the plotted quantities 
remain bounded in agreement with the theory ( Constantin et al (8J.) See figures (7.1a-d). 

(viii) Time periodic attractors. 0.023 ^ v ^ Q.0251 . 

Another window with time-periodic attractors and at least three period-doubling bifurcations was 
also found. Unlike the previous periodic window the period starts with small values at the end of the 
window - the period at v = 0.023 is x = .26 - and reaches the value x = 2.86 at v = 0.0251 after three 
period doublings. The results from this window are summarized in Figures (8.1a-c) for a representa- 
tive member v = .025. 

(lx) Window of Chaotic Trajectories. ? £ v ^ 0.0225 . 

Chaotic solutions are also found in this subwindow with no recognizable patterns to report. The 
ambiguity of the lower bound of this window remains unresolved due to the computational cost of the 
required tuns but also due to the well-known fact that this is not the last chaotic window to be found. 
Window and subwindow lengths become vanishingly smaller thus halting any rational numerical pro- 
gress. The discovery of various attractors becomes a hit or miss problem and, for example, in the 
range .008 £ v £ .0225 we managed to find a strange attractor (v = .0225), a time-periodic attractor 
(v = .017), and bimodal steady attractor (v = .015) before another strange attractor is obtained at 
v = .008. These results are summarized in Figures (9.1a-e). For graphical clarity, the phase-planes in 
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thls set of Figures were obtained by plotting (E , dE/dt ) pairs which are separated by a time which is 
a large multiple of the time-step. 

5. Conclusions. 

We have carried out an extensive numerical study of the Kuramoto-Sivashinsky equation. Our 
emphasis is on the characterization of transition to chaos and we restricted ourselves to odd-parity 
solutions derived from the same initial conditions and which are easily reproducible by a single com- 
puter run. The results for large values of the parameter v, i.e. before chaos or multiple period- 
doubling bifurcations, are completely in line with those of previous investigators (see Introduction). 
The new behavior we found, and which to our knowledge has not been reported before for the 
Kuramoto-Sivashinsky equation, is the computation of several windows which support time-periodic 
orbits and which themselves contain other subwindows which are separated by period-doubling bifurca- 
tions in parameter space. The third periodic window we found (see Section 4 above) provides substan- 
tial numerical evidence of a complete sequence of period-doubling bifurcations before chaos sets in. 
Based on our numerical solutions we computed the Feigenbaum number following the first eight bifur- 
cations of the cascade, and found excellent agreement with Feigenbaum’s universal constant. 
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7. Tables. 

Overview of (lie attracting manifolds 

Window range Description of the attractors 

1 ^ V < oo Constant states. 

.2475 ^ v < 1 Fully modal steady attractors. 

.0756 £ v £ .2472 Bimodal steady attractors. 

.06697 £ v £ .0755 Fully modal steady attractors. 

.0599 £ v £ .06695 Trimodal steady attractors. 

.055238 £ v £ .05985 Periodic attractors. 

.03965 £ v £ .055235 Fully modal steady attractors. 

.03735 £ v £ .0396 Periodic attractors. 

.0344 <, v £ .037348 Tetramodal steady attractors. 

.029756 5 v $ .0343 Periodic attractors. Com- 

plete period-doubling 
sequence towards v=.029756. 
.0252 £ v <, .029755 Chaotic oscillations. 

.024 <. v £ .0251 Periodic attractors. Com- 

plete period-doubling 
sequence towards v=.02512. 

? £ v £ 0.023 Chaotic oscillations. 


Tabic 1 
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First Periodic Window .0598 £ v <; .055238 

Viscosity v 

Period x 

Figures 

.0598 

1.018 


.0595 

1.03 


.059 

1.05 


.058 

1.08 


.057 

1.13 

(2.1a, b) 

.056 

1.18 


.0559 

1.183 

(2.2a, b) 

.05585 

2.374 


.0555 

2.48 

(2.3a, b) 

.0553 

2.59 


.05525 

2.663 


.05524 

2.729 

♦ 



Table 2 
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Second Periodic Window .03735 £ v £ .0396 

Viscosity v 

Period x 

Figures 

.0396 

0.515 


.03955 

0.533 


.0395 

0.53 


.0393 

0.55 


.039 

0.59 


.038 

0.68 

(4.1a-d) 

.03798 

0.684 


.03795 

1.37 

(4.2a, b) 

.0379 

1J8 

(4.3a) 

.0378 

1.40 


.0377 

1.43 

(4Jb-d) 

.0376 

2.87 

(4.4a-d) 

.0375 

2.95 


.0374 

3.00 


.03738 

3.00 


.03736 

3.01 

(4.5a, b) 

.037355 

1.51 


.03735 

1.51 

(4.5c) 


Table 3 
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Third Periodic Window .0343 £ v £ .0299756 

Viscosity v 

Period t 

Figures 

.0343 

.368 


.0342 

.37 

(6.1a) 

.032 

.44 


.0304 

.89 


.0303 

.89 

(6.1b,c) 

.03007 

.88 

(6. Id) 

.030065 

1.76 


.03006 

1.76 

(6.2a) 

.030015 

1.76 

(6.2b) 

.03 

1.76 

(6.2c) 

.029998 

1.76 

(6.2d) 

.0299975 

3.52 


.029995 

3.52 

(6.3a) 

.02999 

3.52 

(6.3b) 

.029985 

3.52 


.029983 

3.52 


.0299815 

3.52 

(6.3c) 

.029981 

7.04 


.02998 

7.04 

(6.4a) 

.029979 

7.04 


.029978 

7.04 

(6.4b) 

.029976 

7.04 


.0299758 

7.04 


.0299756 

7.04 

(6.4c) 

.029975 

14.08 

(6.5a) 


Table 4 
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36-mode Galerkln Expansion 

Subwindow boundary 

Length 

Ratio 

.034615 



.0303175 

.0042975 


.0300494 

.0002681 

16.0295 

.02998665 

.00006275 

4.27 

.029973005 

1.3645X10" 5 

4.6 

.02997007 

2.935X10 -6 

4.65 

.029969435 

6.35xlO -7 

4.62 

.029969305 

1.3X10" 7 

4.88 


Table 5 
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8. Figures 

Figure (1.1) : Steady attracting modal profiles : 

(a) v = .6 .7 .8 .9 . 

(b) v~=~ .26 .28 .3 .35 .4 .5 . 

Figure (1.2) : Steady attracting profiles : 

(a) v = .2 .22 .24 .248 .28 . 

(b) v = .08 .09 .12 .15 .2 . 

Figure (1.3) : Steady attracting trimodal profiles : 

(a) v = .067 .07 .072 .074 .0755 .08 . 

(b) V = .06 .061 .063 .065 .06695 . 

Figure (2.1) : First periodic window, .0598 £ v £.055238 : 

(a) Phase plane of Energy - v = .057 . 

(b) Energy time series - v = .057 . 

Figure (2.2) : First time-periodic window, .0598 £ v £.055238 : 

(a) Phase plane of Energy - v = .0559 . 

(b) Energy time series - v = .0559 . 

Figure (2.3) : First time-periodic window, .0598 £ v £.055238 : 

(a) Phase plane of Energy - v = .0555 . 

(b) Energy time series - v = .0555 . 

Figure (3.1) : Window (iii), fully modal steady attractors, .06697 £ v £ .0755 

(a) Energy against time - v = .0397 . 

(b) Supremum of u(x,t) against time - v = .0397 . 

(c) Supremum of d/dx of u(x,t) against time - v = .0397 . 

Figure (4.1) : Second periodic window, time-period x = .68 : 

(a) Energy against time - v = .038 . 

(b) Phase plane of Energy - v = .038 . 

(c) First Fourier coefficient against time - v = .038 . 

(d) Tenth Fourier coefficient against time - v = .038 . 

Figure (4.2) : Second periodic window, time-period x = 1.37 : 

(a) Phase plane of Energy - v = .03795 . 

(b) Energy time series - v = .03795 . 

Figure (4.3) : Second periodic window, time-period x = 1.37 to 1.43 : 

(a) Phase plane of Energy - v = .0379 . 

(b) Phase plane of Energy - v = .0377 . 

(c) Phase plane of the first Fourier coefficient - v = .0377. 

(d) Phase plane of the sixth Fourier coefficient - v = .0377. 
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Figure (4.4) : Second periodic window, lime-period % = 2.87 : 

(a) Phase plane of Energy - v = .0376 . 

(b) Time series of sixih Fourier coefficient - v = .0376 . 

(c) Phase plane of the sixth Fourier coefficient - v = .0376. 

Figure (4.5) : Second periodic window, period-halving from t = 3.01 lo 1.51 : 

(a) Energy time series - V = .03736 . 

Phase plane of Energy - v = .03736 . 

Phase plane of Energy - v = .03735 . 

Figure (4.6) : Second periodic window : 

Evolution of the spatial profile in the periodic attractor - v = .039. 

Figure (4.7) : Second periodic window : 

Evolution of the spatial profile in the periodic attractor - v = .0375. 

Figure (5.1) : Window (v), steady tetramodal attractors, .0344 £ v < .037348 : 

(a) Energy time series - v = .0373 . 

(b) Steady tetramodal attracting profile - v = .0373 . 

Figure (6.1) : Third periodic window, first period-doubling : 

(a) Phase plane of Energy - v = .0342 . 

(b) Phase plane of Energy - v = .0303 . 

(c) Energy time series - v = .0303 . 

(d) Phase plane of Energy - v = .03007 . 

Figure (6.2) : Third periodic window, second period-doubling : 

(a) Phase plane of Energy - v = .03006 . 

(b) Phase plane of Energy - v = .030015 . 

(c) Phase plane of Energy - v = .03 . 

(d) Phase plane of Energy - v = .029998 . 

Figure (6.3) : Third periodic window, third period-doubling : 

(a) Phase plane of Energy - v = .029995 . 

(b) Phase plane of Energy - v = .02999 . 

(c) Phase plane of Energy - v = .0299815 . 

Figure (6.4) : Third periodic window, fourth period-doubling : 

(a) Phase plane of Energy - v = .02998 . 

(b) Phase plane of Energy - v = .029978 . 

(c) Phase plane of Energy - v = .0299756 . 

Figure (6.5) : Third periodic window, fifth period-doubling : 

(a) Phase plane of Energy - v = .029975 . 

(b) Phase plane of Energy - v = .02997 . 

(c) Phase plane of Energy - v = .02995 . 
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(d) Phase plane of Energy - v = .0299 . 

Figure (7.1) : Chaotic windows, .0252 £ v £ .029755 : 

(a) Energy time series (Chaotic case) - v = .028 . 

(b) Energy time series (Chaotic case) - v = .027 . 

(c) Energy time series (Chaotic case) - v = .0255 . 

(d) Energy time series (Chaotic case) - v = .0252 . 

Figure (8.1) : Time-periodic attractors, .023 £ v £ .0251 : 

(a) Energy time series (Periodic case) - v = .025 . 

(b) Energy time series (Blow-up) - v = .025 . 

(c) Phase plane of Energy - v = .025 . 

Figure (9.1) : Chaotic windows, ? £ v £ .0255 : 

(a) Phase plane of Energy (Chaotic case) - v = .0225 . 

(b) Phase plane of Energy (Periodic case) - v = .017 . 

(c) Energy time series (Periodic case) - v = .017 . 

(d) Energy time series (Bimodal steady attractor) - v = .015 . 

(e) Phase plane of Energy (Chaotic case) - v = .008 . 
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Figure (4.7) 


Evolution of the profile In a periodic attractor 
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Figure (6.2) 
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